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Abstract: A new implementation of the POWHEG method [I] into the Monte-Carlo event 
generator SHERPA [2] is presented, focusing on processes with a simple colour 
structure. Results for a variety of processes, namely e + e _ — !> hadrons, deep- 
inelastic lepton-nucleon scattering, hadroproduction of single vector bosons and 
of vector boson pairs as well as the production of Higgs bosons in gluon fusion 
serve as test cases for the successful realisation. The algorithm is fully automated 
such that for further processes only virtual matrix elements need to be included. 



1 Introduction 

Higher-order QCD corrections by now form an important ingredient to many phenomenological studies and 
experimental analyses at both the Tevatron and the LHC. The impact of these corrections has been similarly 
important for various studies of HERA and LEP data. Calculations invoking such corrections, typically at 
next-to leading order (NLO) in the perturbative expansion in the strong coupling a s , and in very few cases 
also at next-to-next-to leading order (NNLO) accuracy, have been used for a wide range of precision tests of 
our understanding of QCD and the Standard Model. They are also important for the subtraction of back- 
grounds in searches for new physics. When being compared to such calculations, experimental measurements 
are usually corrected for detector effects, while the perturbative result is corrected for hadronisation. Only 
after performing these corrections, theoretical predictions and experimental data are on the same footing. 
Hadronisation corrections are typically determined by using multi-purpose Monte-Carlo event generators, 
such as PYTHIA/PYTHIA 8 [3J, Herwig/Herwig++ gl, or, more recently, SHERPA |2]. In the past two 
decades, such event generators have been the workhorses of particle physics phenomenology, including in 
their simulation many features in addition to the perturbative aspects of a collision, such as hadronisation, 
the underlying event, hadron decays etc.. But, to describe perturbative QCD, they typically rely on leading- 
order (LO) matrix elements only, combined with parton showers, which in turn model QCD radiation effects 
in a leading logarithmic approximation. 

Improvements to this approximation can be obtained through merging methods, pioneered in [5J [B] , and 
further worked out in different varieties at different accuracies and for different parton showers, e.g. in 
[3 15]. In this merging approach, tree-level matrix elements for processes with different jet multiplicities are 
combined with parton showers, avoiding problems related to double counting of emissions. Lately, a new 
formulation has been proposed which can be proved to preserve the formal accuracy of the parton shower, 



independent of the process under consideration Despite varying degrees of formal accuracy amongst 
the various methods, their respective predictions tend to agree on a level expected from such improved 
leading order perturbation theory In most cases the approach leads to a dramatic improvement in the 
description of hard QCD radiation, which makes it a state-of-the-art tool for many analyses, that depend 
on the shape of distributions related to hard QCD radiation. 

However sophisticated at modelling multi-jet topologies, the above methods lack the precision of higher- 
order calculations regarding the overall normalisation, i.e. the cross section of the inclusive process under 
consideration. To accurately estimate uncertainties intrinsic to perturbative calculations, which manifest 
themselves for instance in uncertainties related to choices of renormalisation and factorisation scales, full 
next-to leading order corrections are a sine qua non. 

So far, two different methods to achieve the systematic inclusion of complete NLO corrections for a fixed 
multiplicity have not only been described but also implemented, asserting their practicality. The first one 
has been dubbed Mc@NLO [12]. By now it has been applied to a variety of processes using the framework 
of both the HERWIG [T3] and the PYTHIA [H] event generators, proving its versatility. It effectively relies 
on using the parton shower to perform the subtraction of infrared divergences in the real-emission part of 
NLO corrections, resulting in a separation of the event generation into two parts, one starting from the 
Born kinematics of the hard process, the other starting with real-emission kinematics, i.e. including an extra 
parton in the final state. The weights for these two samples are adjusted accordingly and yield the total 
cross section and the hardest emission, correct to first order in a s . However, this technique may produce 
a number of events with a negative weight. The second method to include NLO corrections into parton 
showers is known as the POWHEG technique pQ. In contrast to Mc@NLO, it can ensure that only events with 
positive weights are generated. To achieve this, the simulation starts with a Born-level event, reweighted to 
include the full NLO correction, i.e. including virtual and real corrections, integrated over the real-emission 
subspace. The first emission is then produced using the exact real-emission matrix element, thus providing 
not only the correct differential cross section, but also the correct radiation pattern to first order in a s . The 
big advantage of this technique is that it can straightforwardly be implemented in a shower-independent 
way, which has been used to some extent in [151 116) . where different processes have been treated in the 
framework of PYTHIA and/or HERWIG. In addition, there have been some implementations in the framework 
of HERWIG H — h [17, 18JHL In fact, because the POWHEG method is very similar to traditional matrix- element 
reweighting [19] , it is relatively simple to implement in a process- independent way and thus very well suited 
for the systematic inclusion of NLO QCD corrections to arbitrary processes. This is reflected by the fact 
that a corresponding toolbox has already been advertised in [5D] . 

Obviously, the ultimate goal is to have a multi-jet merging prescription, accurate at NLO, and to simul- 
taneously respect the logarithmic accuracy of the parton shower. A first step into this direction has been 
presented for the case of e + e _ annihilation into hadrons in [3T]. An alternative approach was suggested 
in [35], extending the method of [5] to NLO accuracy for the core interaction. This second method relies 
on the POWHEG technique to ensure the NLO cross section of the core process. Parallel development has 
brought about a similar algorithm, which will be reported shortly |23] and which is based on the generic 
POWHEG implementation presented here. 

The outline of this paper is as follows: In Sec. [5] the connection between the POWHEG method and parton 
showers will be detailed, introducing also the notation used in this publication. The possibility to implement 
the POWHEG technique on top of an existing parton shower will be discussed. Section [3] briefly introduces 
the SHERPA event generator, which sets the framework for this study. Those parts which are relevant 
for the POWHEG implementation will be presented in more detail. Due to the overall setup of SHERPA, 
incorporating matrix element generators capable of generating the real correction terms for a given process 
and the corresponding differential and integrated subtraction terms in the Catani-Seymour dipole subtraction 
scheme [21], the only missing bit of the NLO calculation driving the POWHEG method are the virtual 
contributions, which have been made accessible through interfaces to BlackHat [25] and MCFM [26,. It is the 
first time that the POWHEG method has been automated using dipole subtraction rather than the Frixione- 
Kunszt-Signer method [27] . The overall quality of the implementation will be exemplified in a number of 
characteristic processes in Sec. 01 including the hitherto unpublished case of W^-pair production in hadronic 
collisions. Sec. [5] contains our conclusions. 

A similar procedure, although restricted to the case of e+e~-annihilation into hadrons was also presented in |10| . 
2 It is worth stressing, that in order to ensure formal accuracy a truncated shower must be used, as pointed out in the 
original publications proposing the method. This option is not available for all parton shower algorithms that have been used 
in the actual implementations. 
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2 Parton showers and the Powheg method 



In this publication, the POWHEG method is reinterpreted as an advanced reweighting technique for standard 
parton showers. The following section introduces the necessary notation and outlines the parallels between 
the POWHEG method and traditional matrix-element reweighting. The starting point of the discussion is 
the factorisation theorem underlying the specific parton-shower model, like the DGLAP equation [55], the 
colour-dipole model [35], Catani-Seymour factorisation [23] or antenna factorisation [3D]. Except in collinear 
factorisation, the splitting functions of the parton shower depend on (at least) one additional parton, which 
is often referred to as the "spectator". In order to make this connection explicit, the notation of a dipole-likc 
factorisation is adopted, which is sufficiently general to discuss all relevant features of the POWHEG method 
and its implementation into the SHERPA event generator. 

2.1 Decomposition of real-emission cross sections 

In the following, sets of n particles in a 2 — > (n — 2) process will summarily be denoted by {a} = {a\, . . . , a n }, 
and the particles will be specified through their flavours{/} = {/i, ...,/„} and momenta {p} = {pi, . . . ,p n }. 
The generic expression for a fully differential Born-level cross section in a scattering process with (n — 2) 
final-state particles can be written as a sum over all contributing flavour combinations as 

da B ({p}) = ^d(T B ({a}) , where da B ({a}) = ^ B ({p}) B({a}) , (2.1) 
{/} 

The individual terms in the sum are given by 

B({a}) =£({<?}) B({a}), B({a}) 

d<f> B ({p}) = ^l^ld0 B ({p}) > £({«};/^ 2 ) 

X\ x 2 

Here, |.M.b| 2 ({a}) denotes the partonic matrix element squared, with all factors due to averaging over 
initial state quantum numbers such as spin or colours absorbed into it, and d<?s({a}) is the corresponding 
differential n-particle partonic phase-space element; S({f }) is the symmetry factor due to identical flavours 
associated to the partonic subprocess, while F({a}) denotes the flux factor and C is the parton luminosity 
given by the corresponding parton distribution functions (PDFs). In the case of leptonic initial states, 
ignoring QED initial state radiation, the parton distribution functions f(x, /i 2 ) are replaced by 5(1 — x). 

In a similar fashion, the real-emission part of the QCD next-to-leading order cross section can be written 
as a sum, this time over parton configurations {a±, . . . , a n _|_i}, i.e. including one additional parton. A 
corresponding subprocess cross section reads 

da R ({a}) = d$ R ({p})R({a}) . (2.3) 



' \M B \ 2 ({S}), 



F{{S}) S({f}) 
x\fh{xi^ 2 ) x 2 ff 2 (x 2 ,H 2 ) ■ 



(2.2) 



At this point, it is helpful to introduce a notation for mappings from real-emission parton configurations to 
Born parton configurations. Such mappings combine the partons a,; and dj into a common "mother" parton 
ay, in the presence of the spectator by defining a new flavour and by redefining the particle momenta. 
To be specific, 

The flavour of the "mother" parton, /y, is thereby fixed unambiguously by the QCD interactions, while the 
flavour of the spectator, remains unaltered. The momentum map guarantees that all partons are kept 
on their mass shell. 

Conversely, any Born parton configuration and a related branching process ij, k — > ij. k determine the parton 
configuration of a real-emission subprocess as 

r^{fu *b\b ; {a}) = [ {f} \lfi\if /j} ■ (2-5) 
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The radiative variables §r\b are thereby employed to turn the n-parton momentum configuration into an 
n+l-parton momentum configuration using the inverse of the phase-space map defined by Eq. (|2.4I) . The 
flavour fj is again determined unambiguously by the QCD interactions. Here, also two obvious generalisa- 
tions of Eq. (|2.4[) shall be defined, &ij,fc({/}) and which act on the parton flavours and on the 
parton momenta only. Correspondingly, such generalisations exist for Eq. (|2.5I) . 

In the soft and collinear limits, the partonic matrix element squared, 1Z({3}), can be decomposed as a sum 
of terms 2?y,fc({o}), 

nm) -> E E ^(H) ■ w 

These terms factorise into a Born-level term and a universal splitting kernel, encoding the transition of to 
a, and a,j |24| . The splitting is associated with a universal procedure for factorising the phase space integral 
into a Born level part and a one-particle radiative phase space. 

d<£ R ({p}) = d$ fl fe({p}))d^j* ({p}) . (2.7) 

The existence of universal decompositions like in Eq. (|2.6[) forms the basis of subtraction methods like the 
Catani-Seymour dipole subtraction antenna subtraction [30 , or the subtraction method of Frixione, 
Kunszt, and Signer [37]. It also serves as starting point for the construction of parton shower algorithms [3T1 
I32j . which aim at approximating parton emissions in the collinear and soft limits of the radiation phase 
space, to resum the associated large logarithms, cf. Sec. 12.21 

However, it is important to stress that, also away from the infrared limits, lZ({a}) can be decomposed into 
a number of terms ~R-ij,k analogous to "D%j,k through 

TZ ijtk ({a}) := PijAWWrn), wh ere P^AM) = ^ -P^(M) (2 8) 

Equation (|2.3[) can now be rewritten as a sum of trivially factorised contributions 
da R ({a}) = E E da B (b tJ , k (M))da%f B ({a}) , (2.9) 



where 

<B««-» = <i«{?})^^y. (2-10) 



and Ry,fc({cf}) — C({a})1Zij t i-({a}). These equations are key ingredients to understanding and implementing 
the POWHEG method. 



2.2 Construction of the parton shower 

Due to the non-Abelian nature of QCD, the terms X^fc in Eq. (|2.6p in general do not factorise on the level 
of squared matrix elements, including all colour contributions. To arrive at a practical model for a parton 
shower, sub-leading colour configurations are therefore neglected, which leads to an assumed factorisation 
on the level of squared matrix elements. In the infrared limits one can then write 

V ijtk {{a}) -+ B(6«,*({3})) Sib f"}Z })) i^-^^K^kipi^Pk) , (2.11) 

Ml/ 1) z PiPj 

where the set of momenta fry, &({??}) is determined by the phase space map of the parton shower modelH 
The quantities )Cij,k are the parton shower evolution kernels, which depend on the parton flavours fi, fj and 
fk and on the radiative phase space. The denominator factor 2pipj or any linearly dependent quantity is 
usually used to define the parton shower evolution variable, in the following denoted by t. 



3 Note that here only parton showers with local energy-momentum conservation are considered. Therefore, the phase-space 
maps {p} R -> {p} B exist. 
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Using the above model, the parton-shower approximation of Eq. (|2.10l) can be derived as 



da 



(PS) ij,k 
R\B 



({a}) = d*«j* ({?}) 



S(h j<k ({f})) 1 



S({f}) ZPiPj 



S7r a s K-ij^k {p% iPjiPk) 



£({a}) 



(2-12) 



Particles produced in the parton shower are resolved at a certain evolution scale and can therefore be 
distinguished from particles at higher and lower scales. At most the two particles <Zj and cij, emerging from 
the same splitting process, can be seen as identical. Hence, the ratio of symmetry factors in Eq. (|2.10[) 
changes to 



g to ({/») 
sttf}) 



i 



1/2 if i, j > 2 and 6, = b } 
1 else 



The integral over the radiative phase space can be written as 
d * l i|s(M) = J^~2 dtdz ^~ J H,k{t-,z,4>) , 



(2.13) 



(2.14) 



with t the evolution variable, z the splitting variable, and <f> an azimuthal angle. Here, J denotes the Jacobian 
factor, that potentially arises due to the transformation of variables. Equation (|2.12l) thus becomes 



da 



(ps) ij,k 

R\B 



At 



a s 1 



({«}) = -rdz — — — Jij,k{t,z,<t>)Kij,k{t,z,<j>) 



2tt 2tt S, 



C(b ijik ({a});t) 



(2.15) 



The assignment of the mother parton, the spectator and the underlying Born process can now be assumed 
to be fixed. Then, the sum runs over all possible real-emission configurations originating from this partic- 
ular Born state instead. Furthermore, assuming independence of the individual emissions, i.e. Poissonian 
statistics, this leads to the constrained no-branching probability of the parton-shower model [33] between 
the two scales t" and t' 



Ag>(^";{a} B )=expJ- £ f ^ P** dz Pgj^(M, 



h=q,g 



1 a. 



x "FT" ^^ij,fe(*)2,i 



£(rrj,k(fi,t,z,<t>', {a});t) 



(2.16) 



£({a};t) j ■ 

It is worth noting that Eq. (|2.16p depends on the underlying Born process, since the flavour and momentum 
of the spectator enter as arguments of J%j k an< i fcij,k- The ratio of C in Eq. (|2.16[) accounts for a potential 
change of the parton luminosity when integrating over the initial-state phase spac^. Note that the partons 
{a}e denote a Born-level set, while in (|2.12p and (I2.15[) {a} denote a set of partons at real-emission level. 
This is also indicated by the subscript B in Eq. (|2.16p . 

Using the definition 



the total cross section in the parton shower approximation reads 



(2.17) 



a * = E/ ^B({p})B({a}) 

+ E E 



{/> 



A^(t ,^;{a}) 



{%,&} ft=9>9 " 

1 a s 
Sij 2ir 



" dt 

to ^ 



dz 



Jij,k(t, Z, (j)) 



• - — — ICij,k(t,z,4>) 



2^ 

^^(/^M^M)^) 

C({a};t) 



(2.18) 



4 Note that, depending on the parton shower model, the Xi do not necessarily fulfil the relation xi = Xi/z |34l 135 
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The scale to acts as the infrared cutoff of the parton shower. Simple inspection shows that the sum in the 
square bracket equals unity, since the second term can be written as 

r' M ^Mim. (2 . 19) 

J t dt 

This makes the probabilistic properties of the parton shower explicit. At the same time it also shows that 
this unitarity leads to the cross section in standard parton-shower Monte Carlos to be exactly the respective 
leading-order cross section. In order to evaluate the formal accuracy of the description of the radiation 
pattern, induced by the second term in the square bracket - the first term encodes the probability that there 
is no resolvable emission off the Born-level configuration - a corresponding observable must be introduced. 
This complicates the discussion somewhat and is therefore postponed to Sees. 12. 41 and 12.51 



2.3 Correcting parton showers with matrix elements 

The aim of this section is to devise a simple method for reinstating 0(a s ) accuracy in the emission pattern 
of the parton shower, i.e. the hardest emission in the parton shower should follow the distribution given by 
the corresponding real-emission matrix element. Loosely speaking, the key idea is to replace the splitting 
kernels K, with the ratio of real-emission and Born-level matrix elements. Thus, instead of the splitting 
kernels, this ratio is exponentiated in the Sudakov form factor and employed in simulating the splitting. 

Comparing Eqs. (|2.10l) and (|2.12[) . a corresponding factor correcting the latter to the former can be easily 
identified. Using Eq. (|2.13l) it reads 

w (m . = d^|B({°» = ?JkPj_ S({/» PlJ , k ({S})n{a}) 

da% S J ij ' k ({3}) 8^ s B{h hk {{a})) /C y -, fe ({a}) ' ^ 

Employing the parton shower approximation, Eq. (I2.1ip . to replace pij t k yields 

"■««}>= >; >; ~ """".'r, :: i m;" 1 ^ — • (2-21) 



v-^ v-* S(b mni i({f})) B(b mni i({a})) 8Tra s 

{m,n}l^m,n UJ J/ 



Note that this implies a corrective weight, which is actually splitter- spectator independent. 

Correcting the parton shower to the full matrix element can thus be achieved through the following algorithm: 

1. Determine an overestimate for Eq. (|2.21D . i.e. find a set of Wy,/ 4 ({3}), such that 

w ( r i] fc(/«' ®R\b\ {a})) < ^/(({J}) for all k and throughout the real-emission phase space. 

2. Replace the parton shower splitting kernels /C^ ^ by W^ j j i ({a}) /Cy- fc. 

3. Accept parton-shower branchings with probability w(r~ <&r\b'i {^}))/W / j5,/ i ({^D- 

It is straightforward to show that the constrained no-branching probability of such a matrix-element corrected 
parton shower reads 

A^ ) (t',t";{a}) = expj- Eifdt H dz f" ^%(t,z,« 

I f™ Jt ' ^ (2.22) 

1 S(r~ ~ k (f i; {f})) R tJ -, fc (r r ^(/ t ,t,z,^;{a})) \ 
X S tj S ({f}) B({a}) J • 

The ratio R/B in Eq. (|2.22p coincides with the ratio in the original publications presenting the POWHEG 
method. In the relatively simple cases treated so far [Tl ll5l[To] . the various symmetry factors in the equation 
above cancel and can be neglected. For more complicated flavour structures this factor may differ from one 
and therefore must be retained. 



6 



Employing again the definition of Eq. (|2.17[) , but this time for the Sudakov form factor constructed from the 
ratio R/B yields the cross section in the matrix element improved parton shower approximation. It reads 



a v = E/ d<MM)B({a}) 
{/} 



A( ME )(t V;{«}) 



E E 

{ij,k} fi=q,g 



16tt 2 



dt 



dz 



2tt 



J ij,k (tf %i 



(2.23) 



1 S{r~ h (f t -{f})) R 4il *(r r7 . £ (/ i ,t,z,0;{^)) 



% 3 s ft V ij,k\ 



S({f}) 



ma}) 



Again, the term in the square bracket equals one and thus reflects the probabilistic nature of this approach. 
Consequently, in the matrix-element improved parton-shower approximation the total cross section is given 
by the Born cross section, although the radiation pattern has improved. For a detailed discussion of the 
real-emission term see Sec. 12.51 



2.4 Approximate NLO cross sections 



In the previous two sections it has become clear that the total cross section of events simulated in a parton- 
shower Monte-Carlo is determined by the "seed" cross section, typically computed at Born level. While 
matrix-element improvement of the naive parton shower picture will lead to radiation patterns which are 
accurate to 0(a s ), the total cross section of the event sample and any observable that can be defined at 
Born level will still be given by the respective leading-order expression. To allow for a simulation with 
next-to-leading order accuracy, including the cross section of the event sample, a prescription to assign a 
corresponding weight and multiplicity of the seed event must be found. 

The solution is to replace the original Born-level matrix element with a modified one pQ, denoted by B, 



da B ({a}) da s ({a}) := d$ B ({p }) B({a}) 



(2.24) 



such that the "seed" cross section, da§, integrates to the full NLO result. When constructing such an 
NLO-weighted differential cross section for the Born configuration, certain approximations must be made, 
since NLO cross sections have two contributions, one with Born-like kinematics and one with real-emission 
like kinematics, both of which exhibit divergent structures. The value of a given infrared and collinear safe 
observable, O, computed at NLO, is given in terms of the Born term B, the real emission term R, and the 
virtual contribution (including the collinear counter-terms), denoted by V, as 



(O) 



(NLO) 



E 

{/} 



d*s({p}) 



B({a})+V({a}) 



0({p}) 



jd$ fi (M)R({a})0(M) 

{/} 



(2.25) 



It is obvious that the real-emission contribution cannot be simply combined with the Born and virtual terms, 
as it depends on different kinematics. In the following, the solution of this problem in the framework of the 
POWHEG method is outlined. 

In order to compute Eq. (|2.25l) in a Monte-Carlo approach, subtraction terms, rendering the real emission 
finite in D = 4 space-time dimensions are introduced. Corresponding integrated subtraction terms regularise 
the infrared divergences of the virtual terms. In the dipole subtraction method [23], the equation above can 
then be written as 



(O) 



(NLO) 



E 



d<MM) 



B({a})+V({a})+I(M) 



0({p}) 



J d$> R ({p}) 

{/} 



R({a})0({p})-J2 E ^AM)0(b t3 , k ({p})) 



(2.26) 
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Note that each Sy,fc defines a separate phase space map and that the observable O in the last term depends on 
bij,k{{p}), rather than {p}, which is a crucial feature of the subtraction procedure. The real and integrated 
subtraction terms Sij^({a}) and I({a}) fulfil the relation 



£ J d$ B ({p})I({a}) = EE E / d*H(W)S y ,*({3}) 

{/} ( l J'} k^ij 



(2.27) 



{/} 



Identifying Dy t with Sj^fc, the term with real-emission kinematics in Eq. (I2.26[) can then be decomposed 
according to Eq. (|2.8|) . resulting in 



0({p}) - S y -, fe ({a}) 0(% fc ({p })) 



{/} {ij} fc #*J' 

In the POWHEG method, this term is approximated as 



(2.28) 



E*«(«)E E 

{/} k ^i-J 



Rii,fc({a}) - S,-j,fe({a}) 



0(ftii,fc({p») 



E d ^(W) E E 



(2.29) 



O(M) 



(r 0i s ({«})) -S«,*(r 5iS ({«})) 

{/} {y,fc}/i=?,9 

This allows the recombination of all contributions to the NLO cross section. Employing Eq. (|2.24p . therefore 
B({a}) = B({a})+V({a}) + !({«}) 



E E /«S 



{ij.fc} A =9,3 



R «,fc( r i5,fe({a})) - s y,fc( r t ~,,fe({a})) 



(2.30) 



In the next section it will be shown that the combination of this term, B with a matrix-element reweightcd 
parton shower yields the attempted 0(a s ) accuracy, not only of the total cross section, but also of the 
real-emission contribution. 



2.5 The POWHEG method and its accuracy 



The key point of the POWHEG method is, to supplement Monte Carlo event samples from matrix-element 
corrected parton showers with a next-to-leading order weight to arrive at full NLO accuracy. This is achieved 
by combining the two methods discussed in Secs. l2.3 land l2.4l To obtain the 0(a s ) approximation to the cross 
section in the POWHEG method, the parton-shower expression of the real-emission probability is combined 
with the approximated initial cross section, dag. This yields the following master formula for the value of 
an infrared and collinear safe observable, O, 



^(POWHEG) = W dMM) B( {5} ) 
{/} 



E E 



i 



16^ 2 j t0 

{y,fe}/i=?,9 



A( ME )(i ,/i 2 ;M)O({p}) 

Jij,k(t, Z, i 



dt 



dz 



o 2vr 

i S '"V, i({/})) R^(r~.,({a})) 



(2.31) 



Si 
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S({f}) 



B({a}) 



A^(t,^;{a})0(r^ k ({p})) 



where obvious arguments of the parton maps r~ g have been suppressed. Clearly, if the observable O on the 
right hand side of Eq. (|2.3ip becomes one, the quantity computed is the total cross section, as for the cases 
discussed in Sees. 12.21 and 12.31 This particular case, however, is insensitive to the details of the radiation 
pattern. To continue the discussion, note that the second term in the square bracket of Eq. (I2.3ip can be 
rewritten as 



E 



dt 



d log Agf(t, 



dt 



A^(t,^ 2 ;{a})0(r lirk ({p})), 



(2.32) 
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which allows to rearrange the expressions in Eq. (j2.31[) as 
(0)( powheg) = d$ B (M)B({a}){o(M) 



{/} 



E 



& ^ A( ME )(^ 2 ;{a}) [0(r~ rk ({p})) - 0({p}) 



Two special cases should now be considered 

• The infrared limit (t —> 0) 

In this case, only the first term in Eq. (|2.33[) contributes, as any infrared safe observable maps the 
real-emission kinematics for collinear (soft) emissions to the kinematics of the (any) underlying Born 
configuration 

o(^i(M))^°o(M) 

The contribution to (£)}( POWHEG ) from this phase-space region is therefore correct to 0(a s ). 

• Hard emissions (t — > ^ 2 ) 

In this case, a cancellation between the first and the last term in Eq. (|2 . 33[) is achieved and only the 
second term remains. Also, A(t, /i 2 ) — > 1 as t — > /i 2 . To 0(a s ) one can then replace B — > B, leading to 



E / d< MM) E E i /W dz / t Ji ^ zA) Sff V ^^(M)) 



Comparing this result with the real-emission term in Eq. (|2.30[) reveals that both expressions differ by 
the factor S({f }) Sij/S(r~ ^({/ })), which arises solely due to the way the real-emission phase space 
is populated by the parton shower (cf. Sec. I2.2[) . 

Therefore, the contribution to ((3}( POWHEG ) from this phase-space region is correct to 0(a e ). 

In the phase-space regions "between" these limits, the POWHEG method interpolates smoothly between the 
two above results. 

3 Realisation of the Powheg method in the Sherpa Monte Carlo 

SHERPA is a multi-purpose Monte-Carlo event generator for collider experiments [5] . The goal of this project 
is a complete simulation of all aspects of the collision. Despite being focused on developments improving the 
treatment of perturbative QCD, over the past years significant improvements have been achieved regarding 
the description of soft QCD and QED dynamics, like the process of hadronisation, the decays of the produced 
hadrons, and the implementation of QED radiation in these decays |36j . One of the traditional key features 
of the SHERPA program, however, is a consistent merging of multi-jet matrix elements at tree- level with the 
subsequent parton shower in the spirit of [5] . An improvement of this method and its consistent implemen- 
tation have been presented in [9] and extended to include hard QED radiation 35 . To this end, SHERPA 
uses its two internal tree- level matrix element generators AMEGIC++ [37] and COMIX [35], which are capa- 
ble of calculating cross sections for processes in the Standard Model (AMEGIC++ and COMIX) and beyond 
(AMEGIC++), involving final states with high multiplicities. Soft and collinear parton radiation is generated 
in SHERPA by means of a parton shower based on Catani-Seymour dipole factorisation [31]. The program 
also allows to steer external modules for the computation of virtual corrections using a standardised inter- 
face [39J. The corresponding real corrections and subtraction terms in the Catani-Seymour formalism [24] 
are then provided automatically by AMEGIC++ [417] , 

SHERPA is therefore perfectly suited to implement the POWHEG method as all prerequisites outlined in Sec. [5] 
are found within a single, coherent framework. In this section the basic features of the corresponding modules 
of the generator are reviewed, as far as they are important in the context of this work. A complete overview 
of SHERPA can be found in . 
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(1) FF 



(2) FI 





a / /Pa 



(3) IF 



(4) II 





Figure 1: Effective diagram for the splitting of (1) a final-state parton connected to a final-state specta- 
tor, (2) a final-state parton connected to an initial-state spectator, (3) an initial-state parton 
connected to a final-state spectator and (4) an initial-state parton connected to an initial-state 
spectator in the standard Catani-Seymour notation. The blob denotes the colour correlated 
leading order matrix element, and the incoming and outgoing lines label the initial-state and 
final-state partons participating in the splitting. 



3.1 Matrix elements and subtraction terms 



For this study, the matrix-element generator AMEGIC++ [37] , is employed. It is based on the construction of 
Feynman diagrams, which are evaluated using the helicity methods introduced in |41j . For the computation 
of NLO cross sections in QCD-associated processes, AMEGIC++ provides the fully automated generation of 
dipole subtraction terms [40] . implementing the Catani-Seymour formalism [24] . As outlined in Sec. 12.41 
such a subtraction procedure is a necessary ingredient to be able and compute NLO QCD cross sections 
with Monte Carlo methods. 

In the Catani-Seymour method, the soft and collinear singularities of the real emission amplitude squared, 
TZ({a}), are removed by a local subtraction term (cf. Eq. (|2.26l) ) 

= EE "WW) 

= E E ^(M)+EE^(H)+EE^ J (M)+EE^(H) ■ 



(3.1) 



i,3 



1,3 a 



a,3 k^j 



a,j b^a 



On the first line, the notation of Sec. l2.4l is adopted, while on the second line the definitions of [24] are restored 
by defining k > 2, a, b = 1,2 and requiring all indices to be mutually distinct. Along those lines, 2?y,fc, 
Ufj, T)°^ and T) a ^ h are the four types of Catani-Seymour dipole terms, as depicted in Fig. [TJ for final-state 
splittings with final-state spectators, final-state splittings with initial-state spectators, initial-state splittings 
with final-state spectators and initial-state splittings with initial-state spectators, respectively. This implies 
that for final state splittings, i.e. i, j > 2, 



S, 



ij,k 



2 ^>ij\* 



and 



-V a . 

2 iJ 



(3.2) 
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Type 


a 


z 




Type 


a 


z 


FF 


Vij,k 


Zl 




IF 


U{ 


•Eik.a 


FI 


1 ^ij,k 


Zi 




II 


Vi 


•Ei,ab 



Table 1: Definition of integration variables in Eq. Q3.5P for the var- 
ious dipole configurations of Fig. [TJ 



while for initial state splittings 

S aj , k = V a J and S aj , b V'^ . (3.3) 

Due to its analytic integrability over the extra emission phase space d$^j^ in D = 4 — 2e dimensions, the 
subtraction term S, in its integrated form 1 of Eq. (|2.27[) , as well as the collinear counter-terms can be added 
back to the virtual contributions to cancel their poles in e, 

{e [V({a}) + l({a})] } ^ = {e [v({a}) + l({d}) + C({a})] } ^ = , (3.4) 

wherein V is the one-loop matrix element convoluted with the Born amplitude and C is the collinear counter- 
term. 

The implementation in SHERPA's matrix element generator AMEGIC++, expanding upon its tree-level capabil- 
ities to generate B and 7Z, is able to generate both the subtraction terms S and their integrated counterparts 
X as well as the collinear counter-term C in an automated fashion. The virtual contributions V, however, 
are obtained from dedicated external codes interfaced using the Binoth-Les Houches Accord [35] ■ Having all 
this at hand, the assembly of the B- function of Eq. (|2.30l) . integrable in D = 4 dimensions, is feasible in an 
automated way. This involves integrating over the real emission subspace of the phase space of the NLO real 
correction to the Born process, cf. Eq. (12.30)) and adding the result to the terms with Born-level kinematics. 

In SHERPA, this integration is performed in a Monte-Carlo fashion, by selecting a single point in the real- 
emission phase space. This technique potentially generates negative weights. In the standard POWHEG 
method, the emergence of such negative weights is suppressed by either performing the integration analyti- 
cally or by sampling over sufficiently many real-emission phase space points. Tests to decide which method 
is better for practical applications are beyond the scope of this publication and will be addressed in a future 
work; here it should suffice to state that, of course, sampling over more than one phase space point is neither 
a conceptual problem nor a practical obstacle. Also, for all processes under study in this publication, no 
problems have been encountered by the possibility of having negative weighted events. Loosely speaking, 
the problem can in no case be more severe than the possibility of having negative weights in a standard NLO 
calculation. Therefore, the only remaining issue is to construct an integration method, which, starting from 
a given Born configuration, is able to fill the real-emission phase space in an efficient manner. Having an 
implementation of the Catani-Seymour subtraction method at hand, the construction of an integrator for 
the real-emission subspace based on CS-subtraction terms is rendered a straightforward exercise. The actual 
integration can be decomposed into three one-dimensional integrals (cf. Eq. (|2.14p ) [2~4] 

^% = 2 -^Aadz^J lhk {a,z^) , (3.5) 

with the two integration variables a and z given in Tab. [TJ cf. |40j . The azimuthal angle, </>, is common to 
all configurations. 

Several different integration channels, each based on a a separate CS dipole, can be combined to yield a 
multi-channel integrator |42j for the real-emission phase space. The a-priori weights in the multi-channel 
can be employed to better adapt to the emission pattern of the process under consideration. Additionally, 
every one-dimensional integrator can be individually improved using the VEGAS algorithm |43) . 
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Type 




Sljfc 


Vij,k 


FF 




Zi 


Vij,k 


FI 


Si 


-Zi 


1 •Eij 1 a 











Type 



IF 



II 



%jk,a 



1 1 - u i 



%jk,a u>j 



i 



Vja.k 



Table 2: Mapping of variables for Eqs. (I3.6[) . (|3.7p and (|3.15[) . Note that the definitions for massless 
partons in [24] are employed. 



3.2 The parton shower 

SHERPA implements a parton shower based on Catani-Seymour (CS) dipole factorisation, which will be 
denoted by CSS [31] . The model was originally proposed in [44] and worked out and implemented in parallel 
in [31] [32]. It relies on the factorisation of real-emission matrix elements in the CS subtraction framework |24) . 
Since the original algorithm has been improved in several ways, which have not yet been compiled in a single 
reference, the parton shower model and its features are briefly reviewed here. 



3.2.1 Ordering parameters and splitting functions 

Consider the process depicted in Fig. [T] where a parton ij, accompanied by a spectator parton k, splits 
into partons i and j, with the recoil absorbed by the spectator k. Conveniently the combined momenta are 
identified as p^ = pi + Pj and Q = pij + pk- A Lorentz-invariant transverse momentum, ky, which acts as 
ordering parameter in the parton shower algorithm, can now be defined as 

£(FS)2 = |q 2 _ m 2 _ m 2 _ m 2| ~__ k z . jk (l _ Z - Jfc ) - (1 - Zjjfe) 2 m\ - z\ ]k m| , (3.6) 

This relation holds, independent of whether the spectator parton is a final- or initial-state particle for all 
final-state splittings, while for initial-state splittings a aj j, again independent of the type of spectator, 
the ordering parameter is given by 

k^ S)2 = \Q 2 - m| - m 2 a - m\\ y ja . k (1 - Zj. ak ) - m| - (1 - z^ ak ) 2 m 2 a , (3.7) 

The precise definition of y and the splitting variables z for the various dipole types are listed in Tab. [5] 

Sudakov form factors for all branching types, taking into account finite masses of final-state partons and 
strictly relying on the Lorentz-invariant variables z and k^, have been derived. The corresponding evolution 
kernels, as defined in Eq. (|2. 1 1|) . read 

ICij, k (^T,z) = ■=- —, — (V ii)fc (k!,z)) , (3.8) 

where A4 P cc is the number of spectator partons of the off-shell particle a™ in the large Nq limit. The 
spin-averaged dipole functions (V), taken in four dimensions, depend on the type of emitter and spectator 
parton and are listed in [31j . Their infrared singularities are regularised through the parton shower cutoff, 
to, typically of the order of 1 GeV 2 . The denominator factor Sij avoids double-counting identical final states 
in final-state evolution and is given by 

A potential shortcoming of the original approach in |31] is that certain dipole functions connecting the initial 
and final state may acquire negative values in some non-singular regions of the phase space. This prohibits 
their naive interpretation in terms of splitting probabilities and leaves the corresponding parts of the phase 
space unpopulated. The problem was solved recently by altering the finite parts of the affected splitting 
functions such that they reproduce corresponding full matrix elements [4"5] . 
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3.2.2 Splitting kinematics 

All branchings in the CSS formalism implement exact four-momentum conservation and the particles are 
kept on their mass-shell before and after every evolution step. The phase-space maps from n- to n+1- 
particle final states are exact and cover the whole phase space. They are, however, not unambiguous. One 
of the most prominent criticisms of this particular parton-shower model is the "unphysical" recoil strategy 
originally employed in configurations with initial-state splitter and final-state spectator. This problem was 
addressed in [33] for the case of massless partons and has been implemented for prompt photon production 
and extended to the fully massive case in [35]. The latter publication defined a general algorithm for 
constructing the emission kinematics, independent of the type of splitter and spectator parton, which is 
summarised as follows: 



1. Determine the new momentum of the spectator parton as 



/ Q 2 + m\ — rrVjj 
Pk= [ Pk 2Q 2 Q 



' \(Q 2 ,s l:j ,m 2 ) 
\{Q 2 ,m 2 ;,ml) 



2Q 2 



Q 



(3.10) 



with A denoting the Kallen function A(a, b,c) — (a — b — c) —Abe and 
Sij = jjij,k (Q 2 - ml) + (1 - yi jt k) (mf + m 2 ) . 



(3.11) 



2. Find the light-like helper vectors I and n as 



I = 



Pij - CtijPk 



1 - ay a* 

where ay = sy/^y^ , a k = m 2 k /li 3 ,k and 
1 



jij^k = 2 In = 



(Q 2 



Pk - Ctk Pij 

1 - Ot,-afc 



ml) + sgn (Q 2 - s i:j - ml) \J X(Q 2 ,Sij,m 2 k ) 



(3.12) 



(3.13) 



3. Express the momenta pi and pj in terms of I, n and a transverse component k±_ as 



Pt = Zi P 



2ln +k± ' 



Pi 



m 2 + k 2 , n ^ 



Zi l in J 1 — z 

The parameters Zi and k^ of this decomposition are given by 

.2 / „~... 2m 2 



2 In 



(Q 2 ,s lj7 m 2 k ) 



i,jk 



Uij,k 



\lij,k\ V 1 — Vij,k 



k 2 



(Q 2 



1-3 " L k 
2™2 =2 2 



;) Vij,k Zi (1 - Zi) - (1 - Zi) m i - Zi rrij , 



(3.14) 



(3.15) 



where £ijf. is defined in Tab. [3] Equations (|3.10p and (|3.15l) are valid for all dipole configurations, i.e. 
initial and final-state branchings with the recoil partner being either in the initial or in the final state. 

To generate an emission using the reweighting technique presented in Sec. 12.31 it is important to be able 
and access matrix-element information during the parton-shower evolution, such that Eq. (|2.21[) can be 
implemented in a process-independent manner. SttERPA provides an interface between its tree-level matrix- 
element generators and its parton shower, which allows for all the necessary interactions. Together with 
an implementation of the phase-space maps 6y j fc({a}), that correspond to the inverse of the above splitting 
kinematics, Eq. (|2.21|) can then be realised easily. 



3.3 Automatic identification of Born zeros 

It was noted in [IB] that Eq. (|2.2ip can develop spurious singularities as the matrix element of the under- 
lying Born process may be zero, while the real-emission matrix element is not. Such configurations do not 
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exponentiate, as 72 is not singular when B — > 0. This fact can be employed to formulate a general solution 
to the problem [IB]. One can split 72 into two parts, a singular one, 7Z^ S \ and a regular one, 72^. 

%M = K Z tt , and 72 (r) = 72 — -— , (3.16) 

Z + H ' Z + H ' v y 

where 

Z = -A- , and ff = K 2 CS -A . (3.17) 

Note that i3 max can be determined during the integration of the seed cross section, while i max is given as a 
universal function of the hadronic centre-of-mass energy, depending only on the definition of t in the parton 
shower model. The resolution factor « res then determines the relative splitting between 72 ^ and 72.W: the 
larger K rcs , the larger the fraction 72 W of 72. 

The necessity of such a splitting of the real emission matrix element can be determined on an event- by- event 
basis by comparing the correction factor to the parton shower of Eq. (|2.20[) . w^fc, to a predefined threshold 
w ijk- Thus, regular non-exponentiated 72^ events are only produced if iOy,fe > w fjk- Such a treatment 
ensures that both the exponentiation of the real-emission matrix element is as inclusive as possible and the 
parton shower correction factor does not get too large, rendering event generation too inefficient. 



4 Results 

This section collects results obtained with the implementation of the POWHEG algorithm in the SHERPA 
event generator. We exemplify the performance in a variety of processes which are listed in Sec. 14.11 with 
their cuts and relevant settings. In Sec. 14.21 the internal consistency of the implementation is checked by 
performing scale variations, cross section comparisons with ordinary NLO calculations, and variations of 
internal parameters of POWHEG. Comparisons of results from the new implementation with predictions from 
tree-level matrix-element parton-shower merging (ME+PS) are presented in Sec. 14.31 Finally, comparisons 
with experimental data are made in Sec. 14.41 



4.1 Process listing 

4-1.1 Jet production in e + e~ collisions 

The annihilation of e + e~ into hadrons is studied at LEP Run 1 energies, Eqms = 91.25 GeV. This setup 
allows to validate the algorithms of Sec.[5]in pure final-state QCD evolution, which is the simple-most testing 
ground. The parton shower cut-off scale has been set to k\ min = 1.6 GeV 2 . Even though the improvements 
in this paper are purely related to perturbative physics, the results are presented after hadronisation with 
the Lund model [3] to make comparison to experimental results more meaningful. The ME+PS samples 
have been generated with up to one additional jet in the matrix elements and the phase space slicing 
parameter was set to log(y cut ) = —2.25. For the virtual matrix elements, the code provided by the BlackHat 
collaboration [25] was used. 



4-1-2 Deep-inelastic lepton-nucleon scattering 

Hadronic final states in deep-inelastic lepton-nucleon scattering (DIS) are studied at HERA Run 1 energies, 
Eqms == 300 GeV. Just like e + e~-annihilation into hadrons, this process boasts a wealth of precise experi- 
mental data. From the theoretical perspective, it is invaluable, as it allows to test QCD factorisation in an 
extremely clean environment. The associated scale, given by the virtuality of the exchanged 7*/Z-boson is 
not fixed, but potentially varies by orders of magnitude, which allows to test perturbative QCD predictions 
in various kinematic limits. Our results are presented at the parton level only, as hadronisation corrections 
have little effect on the observables and the focus lies on the potential improvements on the perturbative 
part of the simulation. The CTEQ6.6 [46] parton distribution functions have been employed and the strong 
coupling has been defined accordingly as a s (mz) — 0.118 with NLO running for both the matrix elements 
and the parton shower. The remaining settings correspond to those in [35]. ME+PS samples have been 
generated with up to one additional jet in the matrix element and the phase space slicing parameters were 
set to Ocut = 5 and Sbis =0.6 (cf. @5]). Virtual matrix elements were provided by BlackHat [2"5] . 
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4-1.3 Drell-Yan lepton pair production 



We investigate Drell-Yan lepton pair production at Tevatron Run 2 energies, simulating pp collisions at 
Ecms = 1-96 TeV. The CTEQ6.6 [46] parton distribution functions are employed and the strong coupling 
is defined accordingly as a s (mz) = 0.118 with NLO running for both the matrix elements and the parton 
shower. A cut on the invariant mass of the lepton pair of 66 < m«/GeV < 116 is applied at the matrix- 
element level. For the ME+PS samples matrix elements with up to one additional jet were generated 
and a phase-space slicing cut of Q C ut = 20 GeV was applied. Virtual matrix elements were provided by 
BlackHat 25 . The factorisation and renormalisation scales for the NLO matrix element were chosen as 
A*ii = I^f = m \ u- I n a ^ tree-level matrix elements SHERPA's default scale choice was employed: The matrix 
element is clustered onto a core 2 — ¥ 2 configuration using a fey-type algorithm with recombination into 
on-shell particles. Scales are defined as the lowest invariant mass or negative virtuality in the core process. 
Hadronisation and multiple parton interactions have been disabled to allow for a study at the parton shower 
level. The Z —> 11 decay is corrected for QED next-to-leading order and soft-resummation effects in the 
Ycnnie-Frautschi-Suura (YFS) approach [36;. 

The three reactions listed in Sec. 14.1. Hl4~T31 essentially amount to one and the same process at the parton 
level, as they only differ by crossing of initial and final state legs. Their combination allows to validate 
the implementation of the matrix-element reweighting in Sec. 12.31 for all possible dipole configurations with 
quark splitters. 

4-1-4 W boson production 

Production of W bosons is presented in pp collisions at Eqms = 1-8 TeV. Although in principle similar to the 
Drell-Yan case, this process is of special interest to validate the automatic decomposition of the real-emission 
term into singular and non-singular pieces, as outlined in Sec. 13.31 If not stated otherwise, the parameters for 
this decomposition are set to K rcs = 4 and wf^ k = 100. The CTEQ6.6 [35] parton distribution functions have 
been employed and the strong coupling has been defined accordingly as a s (mz) — 0.118 with NLO running 
for both the matrix elements and the parton shower. A cut on the invariant mass of the lepton-neutrino 
pair of mt v > 10 GeV was applied at the matrix-element level. For the ME+PS samples matrix elements 
with up to one additional jet were used and a phase space slicing cut of Q cut = 20 GeV was applied. Virtual 
matrix elements were provided by BlackHat [25]. The factorisation and renormalisation scales for the NLO 
matrix element were chosen as fi R = u\ — m 2 ^ , v . In all tree-level matrix elements SHERPA's default scale 
choice was employed, cf. Sec. 14.1.31 Hadronisation and multiple parton interactions have been disabled. 
The W — > Iv decay is corrected for QED next-to-leading order and soft-resummation effects in the YFS 
approach [36j . 

4-1.5 Higgs boson production through gluon-gluon fusion 

The production of Higgs bosons through gluon-gluon fusion is simulated for proton-proton collisions at 
-Ecms = 14 TeV. The coupling to gluons is mediated by a top-quark loop and modeled through an effective 
Lagrangian [37]. Again, this process is technically very similar to the Drell-Yan case, but it also allows to 
validate matrix-element corrections to the remaining initial state splitting functions. Next-to-leading order 
corrections are rather large at nominal LHC energies, with a ratio of K m 2 between the NLO and the 
LO result for the total cross section. This fact has spurred tremendous efforts to perform fully differential 
calculations at NNLO [3Hj and several predictions have been presented which merged such fixed-order results 
with resummation at next-to- next-to-leading logarithmic accuracy TP] , as the process is expected to have 
high phenomenological relevance at LHC energies. However, this publication centres on the behaviour of the 
theory at NLO only, as a prediction beyond this level of accuracy is clearly not within the capabilities of 
the POWHEG method. The CTEQ6.6 [35] parton distribution functions have been employed and the strong 
coupling has been defined accordingly as a s (mz) — 0.118 with NLO running for both the matrix elements 
and the parton shower. A cut for the invariant mass of the r pair of 115 < m TT /GcV < 125 was applied 
at the matrix-element level. For the ME+PS merged samples matrix elements with up to one additional 
jet were used and a phase-space slicing cut of Q cut — 20 GeV was applied. The virtual matrix elements 
have been implemented according to [T5] . The factorisation and renormalisation scales for the NLO matrix 
element were chosen as n 2 R = fi R = m\ . In all tree- level matrix elements SHERPA's default scale choice 
was employed, cf. Sec. 14.1.31 Hadronisation and multiple parton interactions have been disabled. The 
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e + e — > hadrons 


e + p -> e + + j + X 


E cms = 91.2 GeV 


E cms = 300 GeV 
Q 2 > 150 GeV 2 


jU = fJLR = \l F 


Factor 


POWHEG 


NLO 


POWHEG 


NLO 




1/2 


30179(18) 


30195(20) 


3906(9) 


3908(10) 




1 


29411(17) 


29416(18) 


4047(10) 


4050(11) 


2 


28680(16) 


28697(18) 


4180(10) 


4188(11) 



Table 3: Cross sections in pb for e + - e~ annihilation into hadrons at LEP 
and deep-inelastic positron-proton scattering at HERA as calcu- 
lated in the POWHEG framework and in a conventional fixed order 
NLO calculation [40] . 



h — )■ tt decay is corrected for QED soft-resummation and approximate next-to-leading order effects in the 
YFS approach [36]. 

4-1.6 Z -pair production 

The production of pairs of Z bosons is studied for proton-proton collisions at -Ecms — 14 TeV. This is an 
important background for the golden-plated Higgs boson discovery mode at the LHC. Detailed studies of 
the decay properties of the Z bosons and their correlations are known to allow for a determination of some 
properties of the Higgs boson, when found. Among these correlations are, e.g. the relative orientation of the 
decay planes of the bosons. 

The CTEQ6.6 [46 parton distribution functions have been employed and the strong coupling has been defined 
accordingly as a s (mz) = 0.118 with NLO running for both the matrix elements and the parton shower. A 
cut on the invariant mass of each lepton pair of 66 < ma /GeV < 116 was applied at the matrix-element 
level. For the ME+PS samples matrix elements with up to one additional jet were used and a phase-space 
slicing cut of Q cu t = 20 GeV was applied. Virtual matrix elements were provided by MCFM [26] [50] . The 
factorisation and renormalisation scales were chosen as fi R = [ip = m zz . Hadronisation and multiple 
parton interactions have been disabled to allow for a study at the parton shower level. Each Z — > £1 decay 
is corrected for QED next-to-leading order and soft-resummation effects in the YFS approach [36]. 

4-1.7 W + W~ -production 

W^ + W~-production is also studied for proton-proton collisions at -Ecms = 14 TeV. It is worth noting that 
this process hitherto has not been treated in the POWHEG approach. Similar to the Z-pair production, it is 
an important background to the search channel for the Standard Model Higgs boson, at masses around and 
above 130 GeV. Again, in order to suppress this background, distributions which depend on correlations of 
decay products of the VF's in phase space are heavily employed. 

In the simulation here, again the CTEQ6.6 |46j parton distribution functions have been employed and the 
strong coupling has been defined accordingly as a s (mz) = 0.118 with NLO running for both the matrix 
elements and the parton shower. A cut on the invariant mass of each lepton-neutrino pair of mg v > 10 GeV 
was applied at the matrix-element level. For the ME+PS samples matrix elements with up to one additional 
jet were used and a phase-space slicing cut of Q cu t = 20 GeV was applied. Virtual matrix elements were 
provided by MCFM [56] [56] . The factorisation and renormalisation scales were chosen as /j, r = ^l 2 f — Tn^ w . 
Hadronisation and multiple parton interactions have been disabled to allow for a study at the parton shower 
level. Each W — > £v decay is corrected for QED next-to-leading order and soft-resummation effects in the 
YFS approach [36]. 
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pp^ W++X 


pp -> Z + X 


pp~> h + X 


E cms = 1.8 TcV 
m& > 10 GeV 


E cms = 1.96 TeV 
66 < m ee < 116 GeV 


E cms = 14 TeV 
115 < m TT < 125 GeV 


M = Mi? = Mf 


Factor 


POWHEG 


NLO 


POWHEG 


NLO 


POWHEG 


NLO 


mtu/mtt 


1/2 


1235.4(5) 


1235.1(1.0) 


243.96(14) 


243.84(16) 


2.3153(13) 


2.3130(13) 


1 


1215.0(5) 


1214.9(9) 


239.70(13) 


239.59(16) 


2.4487(12) 


2.4474(13) 


2 


1201.4(5) 


1202.0(9) 


236.72(13) 


236.77(15) 


2.5811(13) 


2.5786(13) 




1/2 


1231.0(5) 


1230.3(1.0) 


243.00(14) 


243.06(16) 


2.2873(13) 


2.2869(14) 


1 


1211.8(5) 


1211.7(9) 


239.01(13) 


238.96(15) 


2.4255(12) 


2.4231(19) 


2 


1198.8(5) 


1199.3(9) 


236.23(13) 


236.13(14) 


2.5623(13) 


2.5620(14) 



Table 4: Cross sections in pb for inclusive e + v c ] and Z[-^ e + e~] production at the Tevatron and 

h[—> t + t~] production via a top-quark loop at the LHC as calculated in the POWHEG framework 
and in a conventional fixed order NLO calculation [4*0] . 



4.2 Tests of internal consistency 

The aim of this section is to provide consistency checks on the different aspects of the POWHEG implementa- 
tion in SHERPA. At first, total cross sections as obtained from POWHEG are compared with the corresponding 
results from a standard NLO calculation. In this case, the public release SHERPA- 1.2.20 serves as the refer- 
ence, which includes an implementation of [40] . Results for e + - e~ annihilation into hadrons and deep-inelastic 
positron-proton scattering are presented in Tab. [3] Numbers for inclusive Z-boson production with decay to 
an electron-positron pair, for inclusive VF-boson production with decay to an electron-neutrino pair, and for 
Higgs-boson production via a top-quark loop with decay into r are listed in Tab. HI The agreement between 
the POWHEG results and those of the standard integration method typically is within a la range as given by 
the respective Monte- Carlo errors. 

To examine differences between POWHEG and a parton-shower Monte Carlo regarding the exponentiation of 
the real-emission matrix elements in POWHEG, 1Z can be approximated by 7^ PS ' in Eq. (I2.21j) . Performing 
this replacement does not only constitute a mandatory cross-check, whether the parton-shower approximation 
is retained, but it also estimates the size of corrections that can be expected at all when switching to NLO 
accuracy in the event simulation. Apart from the overall normalisation, in processes with no additional 
phase space dependence introduced by the loop matrix element, the emission pattern in POWHEG should be 
identical to the parton-shower result. This is verified in inclusive Z-boson production at Tevatron energies as 
displayed in Fig. [2] For low transverse momentum (low jet resolution) p± <C mf both distributions coincide 
within statistical errors. For large values the emission phase space is severely restricted in the parton-shower 
approach, as t < up m m z and p± < t. Any contribution to this phase-space region must therefore originate 
from configurations where more than one hard parton recoils against the lepton pair. Such configurations 
are suppressed by higher orders of a s , and therefore the emission rate is gravely underestimated by the 
parton shower. As a direct consequence, all deviations are then manifestations of the exponentiation of 
non-logarithmic terms, which can be sizeable in the hard wide-angle emission region. 

The automatic splitting of the real-emission matrix element into singular and regular contributions as pre- 
sented in Sec. 13.31 contains two unphysical parameters: K rcs , which governs the relative sizes of the exponen- 
tiated, singular part Ttw and the non-exponentiated, regular part TZS T ', and wf^ k , which determines when 
the above separation is actually employed. The effect of n res on the central parton shower reweighting factor 
Wij,k is detailed Fig. [3] There, it can be seen that for values of K rcs chosen neither too low, such that the 
maximum of the reweighting factor rises beyond reasonable bounds rendering the reweighting of the parton 
shower inoperable, nor too high, such that parts of leading logarithmic structure of 7Z are not exponentiated, 
event generation with the accuracy aimed at by the POWHEG algorithm is feasible. Hence, the results of the 
Monte-Carlo simulation should be fairly independent of K Ies and w^ k , if varied within a reasonable range. 
Figure [4] displays predictions for transverse momentum spectra in VF-boson production for several values of 
ftrcs- As expected, no significant variations of the emission pattern can be observed. The small differences 

5 See http://www.sherpa-mc.de 
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Figure 2: — > 1 jet resolution in kx clustered jets and transverse momentum of the e + e~-pair 
in Z /7* boson production at the Tevatron. The standard parton shower effected on 
the leading order matrix elements (red) is compared to the POWHEG formulation (blue) 
and to POWHEG with the real emission matrix element 7Z replaced by its parton-shower 
approximation TZ^ PS ^ (green). 



that can be seen when changing the resolution scale K ros are entirely within the logarithmic accuracy of the 
parton-shower approach and therefore also within the logarithmic accuracy of the real-emission contribution 
in POWHEG. Variations in wf* k only have very little influence on physical distributions. 
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jc res dependence of 10 ^ rt jc res dependence of it^ u ,d 




Figure 3: Dependence of the parton shower correction factor my,fc on the Z-ff-splitting parameter 
K Ies for W~ production at the Tevatron. 
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4.3 Comparison with tree-level matrix-element parton-shower merging 

By comparing POWHEG results to a standard parton shower combined with LO matrix elements (LO+PS), 
it can be established whether observables are produced correctly in regions where the soft/collinear approxi- 
mations in the parton shower are equivalent to the R/B ratios in POWHEG. An example is the distribution of 
the jet resolution scale doi, using the longitudinally invariant A^-algorithm in W/Z + jets production. This 
observable amounts to the fc-r-scale where a 1-jet event is clustered into a 0-jet event. Figure [5] shows that 
there is good agreement between the LO+PS and POWHEG results for doi < 50 GeV. For harder emissions 
the LO+PS approach fails due to the restricted phase space, as explained in the previous section. 
In this publication, we regard the POWHEG method as an advanced matrix-element reweighting technique 
for the parton-shower algorithm; the reweighting is supplemented with local isT-factors to implement full 
NLO corrections. It is therefore useful to compare the respective results with matrix-element parton-shower 
merged samples (ME+PS), which are rescaled by a suitably chosen global if -factor. Such samples are 
known to yield approximate NLO radiation patterns by effectively implementing higher-order matrix-element 
corrections into the parton shower. An implementation of one of the most advanced ME+PS algorithms 
to date is available within the SHERPA framework [S] such that a direct comparison with POWHEG is a 
straightforward exercise. However, because of the lack of virtual contributions in the LO+PS and ME+PS 
samples, an agreement on the total rate cannot be expected. Thus, in the comparisons below the following 
global if-factors were employed: 

• K = 1.038 for e + e~ — > hadrons at LEP energies, 

• K — 1.2 for Z/j* and W production at Tevatron energies, 

• K = 1.2 for ZZ production at the LHC (14 TeV), 

• K = 1.34 for W+W- production at the LHC (14 TeV), and 

• K — 2.1 for Higgs production through gluon fusion at the same LHC energies. 

When comparing POWHEG results to ME+PS results including matrix elements up to the 1-jet final state one 
should obtain a very similar radiation pattern. The observed agreement indeed is very good, as expected. 
Figure [3] shows that, for example, the differential one-jet rates in WyZ-boson production agrees on the 
20% level, even for relatively large scales (doi > 50 GeV). The remaining differences can be attributed to 
the differences in the Sudakov form factors: While POWHEG exponentiates R/B, the ME+PS method uses 
standard Sudakov form factors at the logarithmic accuracy of the parton shower. 

Such differences become visible also in the multiplicity distribution of &t jets with p± > 20 GeV in Drell-Yan 
and W production, cf. Fig. [6j The 0-jet and 1-jet rates agree within 10% between POWHEG and ME+PS, 
but for higher multiplicity final states the POWHEG method predicts significantly more jets. Here a ME+PS 
simulation with more jets in the matrix element would lead to better agreement. 

Now focussing on the properties of the leading jet produced in association with a W or Z boson, the 
transverse momentum of the leading jet is shown in Fig. \S\ Here the LO+PS approach fails to describe the 
hard tail of the distribution, again doe to lacking phase space, while the POWHEG and ME+PS approaches 
agree within 20%. The separation in r]-<f> space between this jet and the W/Z boson is displayed in Fig. [7] 
Clear differences in the shape of the distribution comparing the LO+PS approach with both POWHEG and 
ME+PS are found, as expected, since the other hand, parton showers cover only a restricted area of the 
phase space, and, in addition, they do not encode the full final-state correlations described by the matrix 
elements. Results from the POWHEG and ME+PS methods agree very well, with differences below 10% only. 
A similar finding applies to the transverse momentum of the leading jet, which is shown in Fig. [5] Here 
the LO+PS approach fails to describe the hard tail of the distribution, while the POWHEG and ME+PS 
approaches agree within 20%. 

The transverse momentum of the Higgs boson and the transverse momentum of the leading jet displayed 
in Fig. [9] give a similar picture as in vector boson production: All three methods agree very well for low 
transverse momenta. In the high p± region the POWHEG and ME+PS approaches agree within 15%. 

Figure [TU1 shows that minor differences arise between the LO+PS and the POWHEG and ME+PS approaches 
in the pseudorapidity spectrum of the leading jet. This can be understood as a direct consequence of the 
different transverse momentum distributions in the LO+PS method, as harder jets tend to be more central 
than softer ones. The POWHEG and ME+PS approaches agree well in the central rapidity region and show 
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up to 10% difference only in the forward region. The distribution of r)-(j> separation between the two leading 
jets proves again that the POWHEG and ME+PS predictions are very similar, with deviations below the 
5% level. Again, the LO+PS prediction shows a slightly different behaviour, because of the reasons stated 
above. 

Now we turn to look at diboson production at nominal LHC energies of 14 TeV. Figure [TT] (left) shows 
a comparison of the scalar sum Ht of the transverse momenta of jets and leptons in Z-pair production. 
Deviations of up to 50% become visible between the three compared approaches. This is especially true in 
the high-f/r region. It is well understood that the predictions of the LO+PS approach are softer than either 
of the two other approaches, due to the restricted emission phase space. The relatively large differences 
between the ME+PS approach and POWHEG are naively not expected, but might stem from using consistent 
but somewhat oversimplified scale schemes. This surely should be analysed in more detail, in a forthcoming 
publication, where pair production processes, including WH and ZH would be studied. The transverse 
momentum distributions of the individual Z bosons (Figure [11] right) on the other hand agree very well in 
both approaches, while it is again obvious that the LO+PS sample cannot describe the hard region of this 
spectrum. 

In the azimuthal separation of the two Z bosons, see Fig. 1121 a similar feature as in Ht can be found: The 
events are harder in ME+PS than in POWHEG, leading to increased decorrelation of the boson pair. In 
Figure [12] (right) it can be seen that the angle between the boson decay planes is predicted very consistently 
in all three approaches. 

Properties of the leading jet in W + W~ pair production at LHC energies are displayed in Figure Q21 On 
the left hand side one can see the transverse momentum of the leading jet and on the right hand side the 
separation between lepton and leading jet. For both the ME+PS and POWHEG approaches agree well within 
20% and the LO+PS sample severely underestimates the hardness of the first jet due to its phase-space 
restrictions. 

Figure Q3] displays observables related to the two oppositely charged leptons from the two decays. The 
pseudorapidity difference (left) agrees within 20% for all three approaches, while their azimuthal decorrelation 
is significantly lower in the LO+PS sample than in the ME+PS and POWHEG approaches, which agree very 
well. 
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Figure 5: Predictions for the — > 1 jet resolution in fc*r clustered jets in Z/j* (left) and W (right) boson 
production at the Tevatron. 
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Figure 6: Exclusive jet multiplicity for jets with p± > 20 GeV in the kx algorithm in zT/7* (left) and W 
(right) boson production at the Tevatron. 
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Figure 7: Separation between vector boson and leading jet in zT/7* (left) and W (right) boson production 
at the Tevatron. 
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Figure 8: Transverse momentum of the leading jet in Z/j* (left) and W (right) boson production at the 
Tevatron. 
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Figure 9: Transverse momentum of the Higgs boson (left) and leading jet (right) in gg — > h fusion at 
nominal LHC energies. 
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Figure 10: Rapidity of the leading jet (left) and separation of the leading and second-leading jet (right) 
in gg h fusion at nominal LHC energies. 
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Figure 11: Ht (left) and transverse momentum of the individual Z bosons (right) in ZZ production at 
nominal LHC energies. 
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Figure 12: Azimuthal angle between the two Z bosons (left) and angle between the two Z decay planes 
(right) in ZZ production at nominal LHC energies. 
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Figure 13: Transverse momentum of the leading jet (left) and separation of the lepton and leading jet 
(right) in W + W~ production at nominal LHC energies. 
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4.4 Comparison with experimental data 



The remainder of this section is dedicated to a comparison of results from the POWHEG approach with 
experimental data to assert the improved description of data, provided by this method. 

For the reaction e + e~ — > hadrons at LEP1 energies the LO+PS and ME+PS predictions do not show 
significant differences except in extreme regions of phase space. The POWHEG prediction confirms that 
picture. This is largely due to the fact that the parton-shower algorithm, which is employed in SHERPA is 
based on Catani-Seymour subtraction terms and those terms constitute a rather good approximation to the 
real-emission matrix element in the process e + e~ — > qqg. 

In the distribution of the Durham jet resolution at which 3-jet events are clustered into 2-jet events (Fig. [15] 
left) all three approaches agree very well with the measurement over large parts of the phase space. Only 
in the hard emission region ?/23 > 0.05 deviations from the LO+PS result can be seen. It is encouraging, 
although not surprising to note that both POWHEG and ME+PS describe the data better. Good agreement 
of all three approaches with each other and with the measurement is also observed e.g. for the thrust 
distribution (Fig. [T5l right), the total jet broadening (Fig. [TBI left) and the C-Parameter (Fig. [T5] right). 

As was discussed in [45], deep-inelastic scattering processes offer the opportunity to test perturbative QCD 
in a region where the factorisation scale of the inclusive process, Q 2 , is smaller than the scale of the actual 
event, which might be set e.g. by the transverse momentum of a hard jet. As measurements can be performed 
down to very low values of Q 2 , many hard jets must usually be included in the simulation to achieve a good 
description of data throughout the phase space. This method cannot be used in the context of this work, 
as the POWHEG implementation in SHERPA can so far only be employed for the core process e ± g — > e ± g. 
Therefore, results are presented for the high-Q 2 region only and the discussion of the \ow-Q 2 domain is 
postponed to a forthcoming publication [23] , where the merging of POWHEG samples with higher- multiplicity 
matrix elements will be discussed. Figure [T7] shows reasonable agreement between the POWHEG results and 
experimental data in a measurement of the di-jet cross section performed at the HI experiment |51l 152) . 
Deviations from the LO+PS result are apparent, especially at lower values of Q 2 , as the phase space of the 
parton shower is severely restricted by the low factorisation scale. Similar findings apply to the rapidity 
spectra shown in Fig. [IS] 

The probably most discussed observable probing the influence of QCD radiation in hadron-hadron collisions 
is the transverse momentum of the lepton pair in Drell-Yan production, which is displayed in Fig. 1191 Very 
good agreement, compared to a recent measurement, is found for both the POWHEG and ME+PS approaches, 
while the LO+PS method is not able to describe large parts of the spectrum because of the restricted real- 
emission phase space. The rapidity of the Z boson in Fig. [TO] is described very well by all three approaches. 

The situation is very similar in PF-boson production. A comparison of POWHEG predictions with Tevatron 
data from the D0 experiment [53] is shown in Fig. [50] where very good agreement between the Monte-Carlo 
result and the data can be observed. In addition to the direct comparison the uncertainties related to a 
variation of the renormalisation and factorisation scales are also shown. Thereby, two different strategies are 
pursued: While the inner (dark) band shows the uncertainty related to a variation of the scale in the hard 
matrix elements only, the outer (light) band shows the influence of varying the scales also in the parton- 
shower evolution. It is rather obvious that the latter approach yields the larger variations, as it is associated 
with an uncertainty in the choice of the strong coupling for the real-emission subprocess. While this process 
essentially determines the shape of the transverse momentum distribution in Fig. 1201 it enters at tree-level 
accuracy only, thus leading to a rather large scale dependence. 
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Figure 14: Pseudorapidity difference (left) and azimuthal angle (right) between the two oppositely 
charged leptons in W + W~ production at nominal LHC energies. 
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Figure 15: Durham 2 — > 3 jet resolution (left) and thrust distribution (right) compared to data from the 
ALEPH experiment [54]. 
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Figure 16: Total jet broadening (left) and C-Parameter (right) compared to data from the ALEPH 
experiment [54 . 
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Figure 17: The di-jet cross section as a function of Q 2 in bins of Et.i + £r.2 (left), the three-jet cross 
section as a function of Q 2 (right top) and the ratio of the three- over the two-jet rate as 
a function of Q 2 (right bottom) compared to data from the HI experiment [5TT 152) . 
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Figure 18: The di-jet cross section as a function of rf, compared to data from the HI experiment |51j . 
rj' denotes half the rapidity difference of the two leading jets in the Breit frame. 
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Figure 19: Transverse momentum and rapidity of the Z boson in Drell-Yan lepton-pair production at 
the Tevatron compared to data from the D0 experiment [55] . 
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Figure 20: Transverse momentum of the W boson in W+jets production at the Tevatron compared to 
data from the D0 experiment [S3]- Scale variations of the POWHEG prediction by factors of 
1/2 and 2 are displayed for two different scale schemes, /x^ = [ir = mi v (left) and [If — 
/ir — Tn±j v (right). The inner (dark) band displays the variations associated with redefining 
the scales for matrix elements alone, while the outer (light) band also takes variations in the 
running coupling of the parton shower evolution into account. 
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5 Conclusions and outlook 



In this publication the successful implementation of the POWHEG algorithm into the SHERPA framework 
was reported. The program is fully automated, relying on SHERPA's efficient matrix-element generation 
modules, which allow to construct real correction terms for given processes and their Catani-Seymour dipole 
subtraction kernels in both differential and integrated form. It is worth stressing that this is the first time 
that the POWHEG method has been applied simultaneously to various higher-order calculations using Catani- 
Seymour dipole terms for partitioning the real-emission phase space. This implementation makes a number 
of processes, computed at NLO and available in public program libraries accessible for matching with a 
parton shower. Additional processes are easily added in SHERPA, by merely linking the corresponding code 
for the virtual correction terms^|. 

The implementation was validated by a number of systematic checks, including 

• the stability of cross sections, as exhibited in Tables [3] and 0J 

• the radiation patterns, through comparison with a fake POWHEG algorithm, based on shower kernels, 
c.f. Figure HJ 

• the automated detection of Born zeroes and their stable cure, as shown in Figure [3] and 0] 

• merged LO samples, see Figures [SHHl 

• and comparison with a variety of data, in Figures I15H20I 

It also included, for the first time, the treatment of W^-pair production; results for this, along with some 
plots for Z-pair production are displayed in Figures [TTH14I 

In the near future, more processes with one coloured line only, such as WH and ZH associated production, 
will be added to the SHERPA framework and in a future publication we will also discuss the slightly more 
subtle treatment of processes with more than one coloured line. Furthermore, the methods developed in this 
publication have been used to generate merged samples with an inclusive cross section at NLO accuracy, 
according to the method already presented in |22) ; we will report on our parallel development in a second 
publication [23] . This set of NLO-implementations will be made available in the future, in a new and 
extended release of SHERPA. 
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